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A straightforward algorithm for the symbolic computation of higher-order sym- 
metries of nonlinear evolution equations and lattice equations is presented. The 
scaling properties of the evolution or lattice equations are used to determine the 
polynomial form of the higher-order symmetries. The coefficients of the symme- 
try can be found by solving a linear system. The method applies to polynomial 
systems of PDEs of first-order in time and arbitrary order in one space variable. 
Likewise, lattices must be of first order in time but may involve arbitrary shifts 
in the discretized space variable. 

The algorithm is implemented in Mathematica and can be used to test the 
integrability of both nonlinear evolution equations and semi-discrete lattice equa- 
tions. With our Integrability Package, higher-order symmetries are obtained 
for several well-known systems of evolution and lattice equations. For PDEs and 
lattices with parameters, the code allows one to determine the conditions on these 
parameters so that a sequence of higher-order symmetries exist. The existence of 
a sequence of such symmetries is a predictor for integrability. 
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1 Introduction 



A large number of physically important nonlinear models are completely inte- 
grable, which means that they are linearizable via an explicit transformation or 
solvable via the Inverse Scattering Transform. Integrable continuous or discrete 
models arise in key branches of physics including classical, quantum, particle, sta- 
tistical, and plasma physics. Integrable equations also model wave phenomena in 
nonlinear optics and the bio-sciences. Mathematically, nonlinear models involve 
ordinary or partial differential equations (ODEs or PDEs), differential-difference 
equations (DDEs), integral equations, etc. [||. 

*Research supported in part by the NSF under Grant CCR-9625421. 
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Whichever form they come in, completely integrable equations exhibit analytic 
properties reflecting their rich mathematical structure. For instance, completely 
integrable PDEs and DDEs possess infinitely many symmetries and conserved 
quantities (if the model is conservative). Perhaps after a suitable change of vari- 
ables, the equations have the Painleve property, admit Backhand transformations, 
prolongation structures, or can be written in bi-Hamiltonian form ]|]. 

The existence of an infinite hierarchy of symmetries for integrable equations 
can be established by explicitly constructing the recursion operator that connects 
such symmetries. Finding symmetries and recursion operators for nonlinear mod- 
els is a nontrivial task, in particular if attempted with pen and paper. Computer 
algebra systems can greatly assist in the search for higher-order symmetries and 
recursion operators. 

In this paper we present a direct algorithm that allows one to automatically 
compute polynomial higher-order symmetries for polynomial PDEs in 1+1 dimen- 
sion and polynomial DDEs (semi-discrete equations). The systems of DDEs or 
PDEs must be of evolution type, i.e. first order in (continuous) time. The number 
of equations and the order of differentiation (or shift level) in the spatial variable 
are arbitrary. 

We use the dilation invariance of the given system of PDEs or DDEs to deter- 
mine the form of the polynomial generalized symmetry. Upon substitution of the 
form of the symmetry into the defining equation, one has to solve a linear system 
for the unknown constant coefficients of the symmetry. In case the original system 
contains free parameters, the eliminant of that linear system will determine the 
necessary conditions for the parameters, so that the system admits the required 
generalized symmetry. Our algorithm can thus be used as an integrability test for 
classes of equations involving parameters. 

For the PDE case, a slight extension of our algorithm allows one to compute 
higher-order symmetries that explicitly depend on the independent variables x 
and t. However, in such cases, it is necessary to specify the highest degree of the 
independent variables in the generalized symmetry. 

Once the generalized symmetries are explicitly known, it is quite often possi- 
ble to find the recursion operator by inspection 0. If the recursion operator is 
hereditary then the equation will possess infinitely many symmetries. If the opera- 
tor is hereditary and factorizable then the equation has infinitely many conserved 



quantities |8| |ll[] . 



For Lagrangian systems the set of higher-order symmetries can be shown to 
lead to the set of conservation laws [29, 32]. For equations without Lagrangian 
structure there is no universal correspondence between symmetries and conserva- 
tion laws. For Mathematica algorithms and code to compute conservations laws of 



nonlinear PDEs and DDEs we refer to |13| and 1 14 , [16|] , respectively. The relation- 
ship between symmetries and conservation laws, as expressed through Noether's 
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theorem, is beyond the scope of this paper (see [^] for details). 

The paper is organized as follows. In Section 2 we give the definition of 
generalized symmetry. We also give the three steps of the algorithm to compute 
generalized symmetries of nonlinear evolution equations. The Korteweg-de Vries 
and Boussinesq equations are used to illustrate the computations. In Section 3 we 
extend the algorithm to nonlinear DDEs (the semi-discrete case), with the Toda 
lattice as leading example. The examples in Section 4 are meant to show how the 
technique works for complex equations, like the nonlinear Schrodinger equation 
and one of its discretizations (the Ablowitz-Ladik system). We use a class of 
fifth-order KdV-type equations, and the Hirota-Satsuma system to exemplify the 
method for equations with parameters. In Section 5 we give details about the 
use of our symmetry package in Mathematica. We briefly discuss other software 
for higher-order symmetries in Section 6. Conclusions and an outlook for future 
research are given in Section 7. 



2 Symmetries of Partial Differential Equations 

2.1 Definition 

Consider a system of PDEs in the (single) space variable x and time variable t, 

u t = F(u, u x ,u 2x , ...,u mx ), (1) 

where u and F are vector dynamical variables with the same number of com- 
ponents: u = (m,u 2 , ...,u n ),F = (F 1 ,F 2 , ...,F n ) and u mx = §^£. The vector 
function F is assumed to be a polynomial function in u, u x , u mx . There are no 
restrictions on the order of the system or its degree of nonlinearity. If PDEs are 
of second or higher order in t, we assume that they can be recast in the form (|l[). 

A vector function G(x, t, u, u^, u 2x , ...), with G = [G\, G 2 , G n ), is called a 
symmetry of (|l|) if and only if it leaves (|l|) invariant for the replacement u — > u+eG 
within order e. Hence, 

D t (u + eG) =F(u + eG) (2) 

must hold up to order e on the solutions of (Q) . Consequently, G must satisfy the 
linearized equation || || 

D t G = F'(u)[G], (3) 
where F' is the Frechet derivative of F, i.e., 

F'(u)[G] = |-F(u + eG)| £=0 . (4) 
oe 

In (H|) and (|j) we infer that u is replaced by u + eG, and u nx by u nx + eD"G. As 
usual, Dt and D x are total derivatives. 



U. Gdkt&§ and W. Hereman / Symbolic Computation of Symmetries 



4 



Symmetries of the form G(x, t, u) are called point symmetries. If G(x, t, u, u x , u*) 
they are called classical or Lie-Backlund symmetries, and all others symmetries 
involving higher derivatives than the first are called generalized or higher- order 
symmetries 



The examples used in the description of the algorithm, involve one or two 
dependent variables. For simplicity of notation, the components of u will be 
denoted by u,v,... (instead of u\,U2, etc.). 

2.2 Algorithm 

To illustrate our algorithm, we consider the Korteweg-de Vries (KdV) equation 

ut = 6uu x + u 3x (5) 

from soliton theory. This ubiquitous evolution equation is known to have infinitely 
many symmetries JjJ. 

Key to our method is the observation that (||) is invariant under the dilation 
symmetry (or scaling) 

(t, x, u) -» (X~ 3 t, X~ l x, A 2 u), (6) 

where A is an arbitrary parameter. The result of this dimensional analysis can be 
stated as follows: u corresponds to two derivatives with respect to x, for short, 
u ~ Tpr- Similarly, ~ Scaling invariance, which is a special Lie-point 

symmetry, is an intrinsic property of many integrable nonlinear PDEs and DDEs. 

Our algorithm exploits this scaling invariance to find symmetries, which now 
proceeds in three steps. 

Step 1: Determine the weights of variables 

The weight, w, of a variable is by definition equal to the number of derivatives with 
respect to x the variable carries. Weights are rational, and weights of dependent 
variables are nonnegative. We set w(J^) = 1. In view of @, we have w(u) = 2 
and w{§i) = 3. Consequently, w(x) = —1 and w(t) = —3. 

The rank of a monomial is defined as the total weight of the monomial, again 
in terms of derivatives with respect to x. Observe that (||) is an equation of rank 
5, since all the terms (monomials) have the same rank, namely 5. This property 
is called uniformity in rank. 

Conversely, requiring uniformity in rank for (0) allows one to compute the 



weights of the dependent variables. Indeed, with w(-§z) = 1 we have 



d 

w{u) + w (~Qj.) = 2w(u) + 1 = w(u) + 3, 
which yields w(u) = 2, w(^) = 3. Hence, w(t) = —3, which is consistent with (0). 
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Step 2: Construct the form of the symmetry 

As an example, let us compute the form of the symmetry of rank 7. Start by listing 
all powers in u with rank 7 or less: £ = {l,u,u 2 ,u 3 }. Next, for each monomial in 
£, introduce enough ^-derivatives, so that each term exactly has rank 7. Thus, 

d d 3 d 5 d 7 

— (u 3 ) = 3u 2 u x , -q^^ 2 ) = ^ x u 2x + 2uu 3x , -q^e( u ) = U 5x, (1) = 0. 

Then, gather the resulting (non-zero) terms in a set TZ = {u 2 u x ,u x u 2x ,uu 3x ,U5 X }, 
which contains the building blocks of the symmetry. Linear combination of the 
monomials in 1Z with constant coefficients q gives the form of the symmetry: 

G = ci u 2 u x + c 2 u x u 2x + c 3 uu 3x + c 4 u 5x . (7) 



Step 3: Determine the unknown coefficients in the symmetry 

We determine the coefficients q by requiring that @ holds on the solutions of 
([[]). Compute D t G and use ([j]) to remove Ut, Ut x , ^txx, etc. For given F, compute 
the Frechet derivative (||) and, in view of @, equate the resulting expressions. 
Treating the different monomial terms in u and its x-derivatives as independent, 
the linear system for the coefficients q is readily obtained. 

For (||), we perform this computation with F = 6uu x + us x and G in ([?]). 
Considering as independent all products and powers of u,u x ,u xx , in 

(12ci - \Sc2)u 2 x U2 X + (6ci - \Sc^)uu\ x + (6ci - I8cs)uu x us x + (3c 2 - 60c4)-u| x + 
(3c 2 + 3c 3 - 90a)u 2x U4 X + (3c 3 - 30c 4 )u x u 5x = 0, (8) 

we obtain the linear system for the coefficients Cj : 

S = {12ci-18c 2 = 0, 6ci-18c 3 = 0, 3c 2 -60c 4 = 0, 3c 2 +3c 3 -90c 4 = 0, 3c 3 -30c 4 = 0}. 

The solution is & = ^ = ^ = c 4 . Since symmetries can only be determined up 
to a multiplicative constant, we choose c\ = 30, c 2 = 20, c 3 = 10 and c 4 = 1, and 
substitute this into (|7|). Hence, 

G = 30u 2 u x + 20u x u 2x + 10mi 3x + u^, x . 

Note that ut = G is known as the Lax equation, which is the fifth-order PDE in 
the completely integrable KdV hierarchy. 

Analogously, for (|5|) we computed the (x—t independent) symmetries of rank < 11. 
They are: 

= u x , G (2) = 6uu x + u 3x , G (3) = 30u 2 u x + 20^^^ + Wuu 3x + u 5x , 
G (4) = 140u 3 u x + 70n 3 + 280uu x u 2x + 70u 2 u 3x + 70u 2x u 3x + A2u x u Ax 
+14wn 5a; + u 7x , 
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G® 



630u 4 u x + 1260uu 3 x + 2520u 2 u x u 2x + 1302u x v% x + 420u 3 u 3x 
+966u 2 u 3x + 1260uu 2x u 3x + 756uu x ii4 X + 2h2uj, x u^ x 



+126u u 5x + I68u 2x u 5x + 72u x u ex + 18uu 7x + u 9x . 

These results agree with those listed in the literature (see e.g. @, |3Q|). 
Remarks. 



(i) The recursion operator [3C, p. 312] for the KdV equation is given by 

U = D 2 + 4u + 2u x D~ l . 



(9) 



This operator is hereditary 11] and connects the above symmetries. 
For example, 

lZu x = (D 2 + Au + 2u x D~ 1 )u x = 6uu x + U3 X , 

TZ(6uu x + u 3x ) = (D 2 + 4u + 2u x D~ 1 )(§uu x + u 3x ) 



30u u x + 20u x u 2x + 10uu 3x + u 5x , 



(10) 



and so forth. 



(ii) Instead of working with the definition (^) of the symmetry, one could intro- 
duce an evolution equation, 



u r = G(x,t,u,u x ,u 2x , 



(11) 



which defines the flow generated by G and parameterized by the auxiliary 
time variable r. The symmetry can then be computed from the compatibility 
condition of (||) and (|il]): 



D r F(u, u x ,u 2x ,...,u r 



D t G(x,t,u,u x ,u 2a 



(12) 



One then proceeds as follows: As above, determine the form of the symme- 
try G involving the constant coefficients c%. Then, compute D^G and use 
([l|) to remove Ut,u fx , etc. Subsequently, compute D T F and use (11) to re- 
move \i T ,\i TX , etc. Finally, use (^) to determine the linear system for the 
unknown q. Solve the system and substitute the result into the form of G. 

Applied to our example, T> t G is computed with G in (^). Next, (||) is 
used to eliminate all ^derivatives of u from the expression of DfG. Then, 
compute D T F with F in the right hand side of (|5|), and eliminate all r- 
derivatives through ( |Tl"| ) after substitution of (^). Finally, expressing that 
B T F - D t G = leads to (g). 

Although this procedure (see Ito |23]]) circumvents the evaluation of the 
Frechet derivative, it seems more involved than our algorithm which uses 
the definition (0). 
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2.3 Symmetries Explicitly Dependent on x and t 

The KdV equation (H|) has also symmetries which explicitly depend on x and t. 
Our algorithm can be used to find these symmetries provided that we specify the 
maximum degree in x and t. 

As an example, we will compute the symmetry of rank 2 for (|5|), that linearly 
depends on x and/or t. In other words, the highest degree in x or t in the symmetry 
is 1. 

We start with the list of monomials in u, x and t of rank 2 or less: 

C = {1, u, x, xu, t, tu, tv 2 }. 

Then, for each monomial in C, introduce enough x-derivatives, so that each term 
exactly has weight 2. Thus, 

D x (xu) = u+xu x , D x (tu 2 ) = 2tuu x , D 3 x (tu) = tu 3x , D 2 X {1) = D 3 x (x) = D 5 x (t) = 0. 

Gather the non-zero resulting terms in a set 1Z = {u, xu x ,tuu x ,tu3 X }, which con- 
tains the building blocks of the symmetry. Linear combination of the monomials 
in 7Z with constant coefficients q gives the form of the symmetry: 

G = C\ U + C2 XU X + C3 tuu x + C4 tus x . (13) 

Now, determine the coefficients c\ through C4 by requiring that (jH) holds on the 
solutions of (||). After grouping the terms, one gets 

(6ci + 6c 2 - c 3 )uu x + (3c 3 - 18c4)tul x + (3c 2 - c 4 )u 3x + (3c 3 - \Sc4)tu x u 3x = 0, 
which yields 

S = {6ci + 6C2 — C3 = 0, 3c3 — I8C4 = 0, 3c2 — C4 = 0}. 

The solution is ^ = 3c 2 = ^- = 04. We choose c\ = |, c 2 = g, C3 = 6 and C4 = 1, 
and substitute this into (^). Hence, 

2 1 

G = -u + -xu x + 6ira x + tn 3:r . 

Similarly, we computed other symmetries of (||) that linearly depend on x and 
They are of rank and 2 : 

G = 1 + 6tii x , and G = 2u + xu x + tuj = 2u + xu x + 6tuu x + t?/3 X . 



Our results agree with those in the literature [28] 
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2.4 Example: Boussinesq equation 

For scaling invariant systems such as (||), it suffices to consider the dilation sym- 
metry on the space of independent and dependent variables. For systems that 
are inhomogeneous under a suitable scaling symmetry, such as the example given 
below, we use the following trick: We introduce one (or more) auxiliary param- 
eter^) with an appropriate scaling. These extra parameters can be viewed as 
additional dependent variables, however, their derivatives are zero. By extending 
the action of the dilation symmetry to the space of independent and dependent 
variables, including the parameters, we are able to apply our algorithm to a larger 
range of polynomial PDE systems. 
Consider the wave equation, 

u t t - U2x + 3uu 2x + 3ul + au 4x = 0, (14) 

(a constant) which was proposed by Boussinesq to describe surface water waves 
whose horizontal scale is much larger than the depth of the water |Q] . 

To apply our algorithm, we must first rewrite (O) as a first-order system, 



ut = v x , 

vt = u x - 3uu x - au 3x , (15) 

where v is an auxiliary dependent variable. It is easy to verify that the terms 
u x and au 3x in the second equation obstruct uniformity in rank. To circumvent 
the problem we introduce an auxiliary parameter [3 with (unknown) weight, and 



replace (15) by 



u t = v x , 

v t = (3u x - 3uu x - au 3x . (16) 

As described in Step 1 we compute the weights from 

d 

w(u) + w(-^) = w(v) + 1, 
d 

w(v) + w(-) = w(j3) + w(u) + 1 = 2w(u) + 1 = w(u) + 3. 
at 

This yields 

d 

w(u) = 2, w{v) = 3, w{(3) = 2, and w{— ) = -w{t) = 2, 

and the scaling properties of (|l~6|) are u ~ (3 ~ ^ ~ v ~ J^-. Indeed, (|i~6|) is 
invariant under the dilation symmetry 

(x, t, u, v, (3) -» (A _1 x, \~ 2 t, A 2 u, X 3 v, X 2 (3). 



U. Gdkt&§ and W. Hereman / Symbolic Computation of Symmetries 



9 



Observe that all the monomials in the equations in (16) have rank 4 and 5. There- 
fore, for any symmetry G of (|l6|), 

rank(C?2) = rank(Gi) + 1 = rank(Gi) + w(v) — w(u). 

Let us construct the form of the symmetry G = (G±, G2) with rank(Gi) = 6 and 
rank(G2) = 7. First, list all monomials in u,v and f3 of rank 6 (respectively rank 
7) or less: 

L\ = {1, /3, (3 2 , /3 3 , u, f3u, p 2 u, u 2 , /3u 2 , u 3 , v, (3v, uv, v 2 }, 

£2 = {1, /?, /3 2 , /3 3 , u, (3u, p 2 u, u 2 , f3u 2 , u 3 , v , f3v, (3 2 v, uv , f3uv } u 2 v , v 2 }. 

Next, for each monomial in C\ and £2, introduce the necessary x-derivatives, so 
that each term in L\ exactly has rank 6, and each term in C2 has rank 7. Keeping 
in mind that (3 is constant, and proceeding with the rest of the algorithm, we 
obtain: 

(1) 2 
G\ = u x v + uv x + -av 3x , 

2 2 

G$p = f3uu x - 3u 2 u x + vv x - 6au x u 2x + -a(3u 3x - 3auu 3x - -a 2 u 5x . (17) 

Finally, setting (i = 1 in (|17|), one obtains a symmetry of (|l^) although initially 
this system was not uniform in rank. We list one more higher-order symmetry of 

©: 

3 2 8 

G? = u Ux --u 2 u x + vv x -5au x U2 X + -au 3x -2auu 3x --a 2 U5x , 

G? = uv x + vu x -3uu x v- 3 -u 2 v x -2a U 2 X v x -3au x V2 X -au 3x v + 2 -av 3x 
8 2 

-2auv 3x - —a v 5x . 
15 



3 Symmetries of Differential-difference Equations 

3.1 Definition 

Consider a system of DDEs, 

u n = F(...,u n _i,u n ,u n+ i, ...), (18) 

where the equations are continuous in time, and discretized in the (single) space 
variable. As before, u ra and F are vector dynamical variables with any number 
of components, and F is assumed to be a polynomial with constant coefficients. 
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There are no restrictions on the level of shifts or the degree of nonlinearity. If 



DDEs are of second or higher order in t, they must be recast in the form (18). 

A vector function G(..., u n _i, u n , u n+ i, ...) is called a symmetry of (|]) if the 
infinitesimal transformation 



u n — > u n + eG(..., u n _i, u n , u ra+ i, ... 

leaves (|l8| ) invariant within order e. 

Consequently, G must satisfy the linearized equation 

D t G = F / (u n )[G], 

where F' is the Frechet derivative of F, defined as 

F'(u„)[G] = ^F(u n + 6G)| e=0 . 



(19) 



(20) 



(21) 



Of course, (19) means that u n+ k is replaced by u n+ fc+eG| n ^ n+fc . For compactness 
of notation, in (20) and (21) we used F'(u n ) instead of F'(..., u n _i, u n , u n+ i, ...). 



Also for notational simplicity, in the description of the algorithm below, the 
components of u n will be denoted by u n ,v n , etc. We use F±, F2, ... and G\, G2, 
to denote the components of F and G, respectively. 

3.2 Algorithm 



As the leading example, we consider the one-dimensional lattice [17, 38] 

ij n = exp (y n _i - y n ) - exp (y n - y n +\), (22) 



due to Toda. In (|22j), y n is the displacement from equilibrium of the nth particle 
with unit mass under an exponential decaying interaction force between nearest 
neighbors. With the change of variables, 



}Ju 



exp (y n -y n+ i), 



the Toda lattice (22) can be written in polynomial form 

u n = w n _i - v n , i) n = v n (u n - u n+ i). 
Observe that (^) is invariant under the dilation symmetry 

(t,u m v n ) -> (X^t, Xu n , \ 2 v n ), 



(23) 



(24) 



where A is an arbitrary parameter. Thus, u n corresponds to one derivative with 
respect to i, or u„ ~ and, similarly, v n ~ Jp-. Our 3-step algorithm exploits 
this scaling property to find symmetries. 
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Step 1: Determine the weights of variables 

In contrast to the algorithm for PDEs, we have to define the weight, w, of variables 
in terms of the number of derivatives with respect to t, and we set = 1- 

Weights of dependent variables are nonnegative, rational, and independent of n. 



In view of (24), we have w(u n ) = 1, and w(v n ) = 2. 



The rank of a monomial is defined as the total weight of the monomial, again 



in terms of derivatives with respect to t. Observe that in the first equation of (23), 
all the monomials have the same rank, namely 2, and in the second equation, all 
the terms have rank 3. 

Conversely, requiring uniformity in rank for each equation in (|23|) allows one 
to compute the weights of the dependent variables. Indeed, 

w{u n ) + l = w(v n ), w(v n ) + l = w(u n ) + w(v n ), 

yields w(u n ) = 1, w{v n ) = 2, which is consistent with (p4|). 
Step 2: Construct the form of the symmetry 

As an example, we compute the form of the symmetry of rank (3,4), i.e. G\ and 
G<i will have ranks 3 and 4, respectively. Start by listing all monomials in u n and 
v n of ranks 3 and 4, or less: 

C\ — {u n , u n , u n v n , u n , v n } , 

r — J 4 3 2 2 2 \ 

Next, for each monomial in L\ and L<i, introduce the necessary t-derivatives. so 
that each term exactly has rank 3 and 4, respectively. At the same time, use 
to remove all t— derivatives. Doing so, based on C\, we obtain 



u -. d° 



^ t (ul) = 2u„u„ = 2»„»„_! - 2u n v n , ^™-„„ +1 „„ 



d 2 , v d , . . d 

-r^{u n ) = —(u n ) = —{Vn-l ~ v n ) = U n _iV n -l - U n V n -l - U n V n + U n+ iV n . 

dt dt dt 



Gather the resulting terms in a set: 1Z\ = {u^, u ra _it> ra -i, u n u n _i, u n v n , u n+ iv n }. 
Similarly, based on the monomials in £2, we get 



^2 = {ui^l^Vn^i^n^iUnVn-i^nVn-iiVn^Vn-l^l-iiUnVn, 
.2 „, „, „, „,2 



2 2 \ 

n^n+lVn-, U n+ iV n , V n -\V n , V n , fn^n+l /• 
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Linear combination of the monomials in TZ\ and 1Z 2 with constant coefficients Cj 
gives the explicit form of the symmetry: 

G\ = Cl V? n + C 2 U n -lV n -l + C 3 U n V n -l + C4U n V n + C 5 U n+ lV n , 
G 2 = C G U n + C 7 U n _ 1 f; n _l + C 8 U n -\U n V n -\ + C 9 U n V n -l + Cio V n ~ 2 V n -\ + 
Cn + C12 U^Vn + C13 U n U n+ iV n + Cu U n+1 V n + C15 Un-l^n + 

Cl6^n + Ci 7 V n U n+ i. (25) 

Step 3: Determine the unknown coefficients in the symmetry 

To determine the coefficients q we require that (^0|) holds on any solution of 
(lilf). Compute D t G and use (ll8|) to remove all u n -i, u n , u n+ i, etc. Compute 
the Prechet derivative ( pl| ) and, in view of (p0|), equate the resulting expressions. 
Considering as independent all the monomials in u ra and their shifts, we obtain 
the linear system that determines the coefficients c%. 

Applied to ( |23| ) with (|25|), we obtain the solution 

ci = c 6 = c 7 = c 8 = c 9 = cio = cu = ci 3 = ci6 = 0, (26) 

-c 2 = -c 3 = C 4 = C 5 = -C12 = C14 = -C15 = C17. (27) 

Therefore, with the choice cyj = 1, the symmetry is 

G\ = U n V n - U n -iV n -i + U n+ iV n - UnV n -l, 

G 2 = u 2 n+1 v n - u\v n + v n v n+ i - v n -iv n . (28) 

It is easy to produce new completely integrable DDEs based on these symmetries. 
For instance, the DDE system 

U n = Gl=U n V n -U n -lV n -l+U n+ lV n -U n V n -l, 

v n = G 2 = u 2 n+1 v n - u 2 n v n + V n V n+ l - V n -lV n . (29) 

is also completely integrable. 

To illustrate the effectiveness of our algorithm to filter out integrable cases 
among systems of DDEs with parameters, consider a parameterized version of the 
Toda lattice, 

u n = a u n _i - v n , v n = v n (/? u n - u n+ i), (30) 

where a and are nonzero constants. In [^] it was shown that (|30| ) is completely 
integrable if and only if a = /3 = 1. 

Using our algorithm, one can easily compute the compatibility conditions for a 
and (3, so that (|3^) admits a polynomial symmetry, say, of rank (3,4). The steps 
are the same as for (p3|). However, the linear system for the Cj is parameterized 
by a and (3 and must be analyzed carefully. This analysis leads to the condition 
a = P = 1. 
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For a = j3 = 1, ( |30| ) coincides with (|23"|), for which we computed symmetries 
with ranks (4,5) and (5,6). They are: 

Gi = u 2 n v n + u n u n+ iv n + u 2 n+l v n + vl + v n v n+ i - Un_ x v n -i - U n -lU n V n -l 
G$p = u n+ ivl + 2u n+1 v n v n+1 + u n+2 v n Vn+i ~ u z n v n + ul +1 v n 

-Un-lVn-lVn ~ 2u n V n -\V n - U n V 2 , 

Gf = u^u n + u 2 n u n+ iv n + u n u 2 n+1 v n + ?4 + it>„ + 2« n ^ + 2n„ + i^ 
-n„v n _ 2 Wn-i - 2n„_i^_ 1 - 2u n v£_ 1 , 

G^ 2) = U* +1 U n - Un v n - U 2 n _ x V n -lV n - 2u n _ 1 U n V n -iV n - 3u 2 n V n -iV n 

-v n - 2 v n -iVn - v 2 n _ x v n - 2u 2 vl + 2*4+1^ - v n -\v 2 n + 3n^ +1 w n f„ + i 

+2-U n+ iU n+ 2f„U„ + l + U 2 n+2 V n V n+ i + V 2 V n+ i + U n «n+1 + v nV n +\V n +2- 



4 More Examples 

4-1 Nonlinear Schrddinger Equation 

The nonlinear Schrodinger (NLS) equation pL 

- fee + 2|<?| 2 <? = 0, (31) 

arises as an asymptotic limit of a slowly varying dispersive wave envelope in a 
nonlinear medium, and as such has significant applications in nonlinear optics, 
water waves, and plasma physics. Together with the ubiquitous KdV equation 
(||), the completely integrable NLS equation is one of the most studied soliton 
equations. 

In order to compute the symmetries of (31) we consider q and q* as indepen- 
dent variables and add the complex conjugate equation to (|3l|). Absorbing i in 
the scale of t, we get 

qt ~ + 2q 2 q* = 0, q* t + q* 2x - 2q* 2 q = 0. (32) 

Since w(q) = w(q*), we obtain 

d 

w(q) = w(q*) = 1, and w (~gj.) = ~ w (t) = 2- 
Hence, the following are the symmetries of ranks (4,4), (5,5), and (6,6): 
Gi = Sqqxq* + q3x, G%> = Gl {1) = -6qq*q* + q% x , 
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Gf = -6q 3 q* 2 + Gq 2 x q* + 4qq x q* + 8qq 2x q* + 2q 2 q* 2x - q ix , G { 2 2) = G\ {2) 
Gf ) = 30q 2 q x q* 2 - 10q x 2 q* - 20q x q 2x q* - Wqq2 X q* x - Wqq x q 2x 
-10qq 3x q* + q 5x , Gf = Gf\ 

4-2 Fifth-Order Korteweg-de Vries Equations 

Consider the parameterized family of fifth-order equations, 

u t + au 2 u x + (3u x u 2x + 7uu 3x + u 5x = 0, (33) 



where a, f3, 7 are nonzero constants. Integrable cases of fl33| ) are well known in the 
literature |9|, 26, |34"f. Indeed, for a = 30, (3 = 20, 7 = 10, equation ( |33[ ) reduces 



to the Lax equation P7|| . The SK equation, due to Sawada and Kotera [35|, and 



Dodd and Gibbon ||, is obtained for a = 5, j3 = 5, 7 = 5. The KK equation, due 
to Kaup [25| and Kupershmidt, corresponds to a = 20,(3 = 25,7 = 10- 

The scaling properties of (]33|) are such that u ~ j ~§i ~ "§x^ ■ Using our 
algorithm, one easily computes the compatibility conditions for the parameters 
a, P and 7, so that (|33| ) admits a symmetry of fixed rank. The results are: 



Rank 5: G = uu x + ^-u 3x is a symmetry of (|33|) provided that 



Rank 3: G = u x . is a symmetry of ( p3| ) without any conditions on the parameters. 

t> = ^7 2 , and/J = 2 7 . (34) 



The Lax equation satisfies (|34|). Since the KdV equation (|5|) is a member of Lax 
hierarchy, condition (^J) comes as no surprise. 

Rank 7: Equation (^) is of rank 7. The stationary part of ( |33[ ) is the symmetry. 
Rank 9: Three branches emerge: 
(i) If condition ( |34[) holds then 



« 5 „ 20 5 2 50 30 

G = u°n a . + H uu x u 2x + -u u 3x + -^u^ar + 



10 50 

H ~MW5a: + 7T^ u 7i 

/yZ I /yO 



(ii) If 

a. = 

5 



a = -7 2 , and /3 = 7 (35) 



holds, one has the symmetry 

3 15 o 45 15 2 225 75 

G = U ttj; + — U x + —UU x U2x + —U U 3x + — -kU2xU3x + ~Z~^U X U^ X 

47 27 27 47^ 2^y z 

75 375 

+ + ^3 U 7x- 
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The SK equation satisfies the condition (|3 

(iii) One has the symmetry 

o 75 3 135 15 2 225 525 

G = U U x + — U x + ——UU x U 2x + TT U M 3x + 7; n u 2x u 3x + ^-^asUfc 
87 47 27 27^ 87^ 



75 375 
+^uu 5x + ^u 7x 

provided that 



1 2 , - 5 



a = -7 , and /3 = -7, (36) 

which holds for the KK case. 

Rank 11: One obtains the symmetry 

4 20 o 40 2 620 2 20 o 460 2 

G = u u x H u< H u + -^u x u 2x + — «°u 3x + — -^^3^ 

7 7 37 z 37 37^ 

200 120 400 20 2 800 

H y-UU 2x U3 x H ^-UU x U ix H o-U3a;W4a: H n^* ^50; + T^T^a^a; 

/y« ry^ ,-yO ,-y.i 



800 200 1000 



UxU6x + 7T^ uu 7x + 7^w«9. 



7 7 3— ox , ?7 3-^ ' 637 4 



a; 



provided that the condition (|34| ) for the Lax hierarchy is satisfied. 

In summary, our algorithm allows one to filter out all the integrable cases 
in the class (33). Alternatively, in [13| we investigated the conditions on the 
parameters a, /?, 7 such that (|33|) admits an infinite sequence (perhaps with gaps) 
of polynomial conservation laws. The conditions in [13] are exactly the same as 
the ones above. 



4-3 Hirota and Satsuma System 

Hirota and Satsuma [^] proposed a coupled system of KdV equations, 

u t = 6auu x + 6vv x + ctU3 X , 

v t = 3uv x + v 3x , (37) 

where a is a nonzero parameter. System ( |37|) describes the interaction of two long 
waves with different dispersion relations. It is known to be completely integrable 
provided a = — |. The scaling properties of j37| ) are such that u ~ v ~ J^, ~ 



J^s. System (|37|) is of rank 5. In a search for the symmetry of higher ranks, we 
obtained the symmetry of rank 7: 

1 2 2 2 1 2 2 2 

&2 = -o"" % - o V V X ~ oU 2x V x ~ ^U x V 2x ~ ~UV Sx - —V 5x , 

3 3 3 3 3 15 
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provided that a = — ^, which is the condition for complete integrability of (|37|). 
4-4 Volterra Chain 

Consider the following integrable discretization of the KdV equation: 

<=«nK+l-«n-l), (38) 

which is also known as the Kac-Van Moerbeke equation. It arises in the study of 
Langmuir oscillations in plasmas and in population dynamics |l], 

Note that ( |38| ) is invariant under the dilation symmetry (t, u n ) — > (\~ 1 t, \u n ). 
Hence, u n corresponds to one derivative with respect to i, i.e. u n ~ 

We computed the symmetries of (|38|) with ranks 3 through 5. They are: 

= U n U n+1 (u n + U n+1 + U n+2 ) - U n _iU n (u n _2 + Un-1 + U n ), 
/q\ q 2 2 3 2 2 

G { > = u n u n+ i + 2u n u n+1 + u n u n+1 + U n U n+ lU n+ 2 + 2u n u n+1 u n+ 2 

+U n U n+ iU 2 n+2 + U n U n+1 U n+2 U n+3 - U n _ 3 U n - 2 Un-lU n ~ U^_ 2 U n -lU n 

2 3 2 3 

— 2u n -2U n _iU n — U n _iU n — U n -2U n -\U n — U n -iU n , 

/q\ ^ 22 32 23 4 3 

G ( ' = u n u n+ i + u n -iu n u n+1 + 3u n u n+1 + 3n n -u n+1 + u n u n+1 + -u n u n+ i-u n+2 
+Au 2 n u 2 l+l u n+2 + 3u n u\ +l u n+2 + ?4ti„ +1 ?4 +2 + 3u 

+-u n n n+ iu n+2 n^ +3 + n n -u n+ i-u n+2 n n+3 -u n+4 - ■u n _ 4 u n _ 3 u n _ 2 -u n _iu„, 
-u 2 _ 3 u n - 2 u n -iu n - 2-u n _3^_ 2 n n _iu n - v? n _ 2 u n -iu n 
—2u n ^ 3 u n — 2 u 2 n __^u n — 3u 2 l _ 2 u 2 l _iU n — 3u n — 2 uf l _iU n — u\_^u n 
-u n -3U n - 2 u n -iu 2 n - u n _ 2 u 2 n _ x u 2 n - \u n - 2 u 2 n _ x u 2 n - 3u\_ x u 2 n 

3 o 2 3 4 2 2 

~ u n—2U n -lU n — 6U n _iU n — U n -\U n — U n _-yU n U Tl J r \. 

Ignoring a trivial misprint in []28|| , Mikhailov et al. listed the symmetry 

4-5 Ablowitz-Ladik Discretization of the Nonlinear Schrodinger Equation 

In |2], [J, Ablowitz and Ladik studied some of the properties of the following 
integrable discretization of the NLS equation: 

iu n = u n+l - 2u n + -u n _i ± u* n u n {u n+1 + u n -i), (39) 

where u* n is the complex conjugate of u n . We continue with the plus sign; the 
other case would be treated analogously. Instead of splitting u n into its real and 
imaginary parts, we treat u n and v n = u* n as independent variables and augment 
( |39"| ) with its complex conjugate equation. Absorbing i in the scale on i, we get 

u n = u n+1 -2u n + u n -i + u n v n (u n+ x + u n -i), 

V n = ~{Vn+l ~ 2v n + V n _i) - U n V n (v n+1 + V n -i). (40) 
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Since v n = u* , we have w(v n ) = w(u n ). Neither of the equations in (|4C|) is uniform 
in rank. To circumvent this problem we introduce an auxiliary parameter a with 
weight, and replace (^0|) by 

u n = a(u n+ i - 2u n + u n _i) + u n v n (u n+ i + it n _i), 

= -a(f n +l - 2w„ + U n _l) - U n V n (v n+ i + U n -l). (41) 

Uniformity in rank requires that 

w(u n ) + 1 = w(a) + w(u n ) = 2w{u n ) + w(v n ) = 3w(u n ), 
w(v n ) + 1 = w(a) + w(v n ) = 2w{v n ) + w(u n ) = 3w(v n ), 

which yields w(u n ) = w(v n ) = ^, w(a) = 1, or, u n ~ v\ ~ a ~ -jr. 

Recall that the 'uniformity in rank' requirement is essential for the first two 
steps of the algorithm. However, after step 2, we may set a = 1. The computations 
now proceed as in the previous examples. We searched for symmetries of ([!(]) of 
ranks (2,2) through (7/2,7/2), and found symmetries of ranks (5/2,5/2) and 
(7/2,7/2). To save space, we only list the symmetries of rank (5/2,5/2) : 

Gi = -u n+2 - u n u n+ iv n ^i - u 2 n+1 v n - u n u n+2 v n - u 2 n u n+ iv n -iv n 

-U n U n+1 V n - U n+ lU n+2 V n+ l - U n U n+1 U n+2 V n V n+ l, 
(1) 2 

G\ = v n - 2 + u n -iv n _ 2 v n -i + u n v n _ x + U n V n ^ 2 V n + U n+ lV n _iV n 
+U n ^iU n V n - 2 V n -iV n + U n V 2 n _ x V n + U n U n+ lV n ^lV n , 

G^ = -U n _2 - U n _ 2 U n -lV n -l - U^-lVn - U n - 2 U n V n - U n - 2 U n -\U n V n -\V n 
-v? n _ x u n v\ - U n -lU n V n+ i - U n -lU 2 n V n V n+ l, 

G^ = U n -lV n V n+ l + U n -lU n vlv n+ l + U n vl +1 + U^VnV 2 ^ 

+V n +2 + U n V n V n+2 + U n+ lV n+ lV n+2 + U n U n+ lV n V n+ lV n+2 . 

4-6 Generalized Toda lattices 

In ]3q| , the integrability of the chain 

y n = y n+ ie (yn+1 - yn) - e 2 (y-+i-^) _ y n _ ie to»-v»-ti + ^(jm-jm-^ ( 42 ) 

which is related to the relativistic Toda lattice has been studied. With the change 
of variables, u n = y n , v n = exp {y n +i — Un), lattice (^) can be written as 

u n = v n (u n+ i - v n ) - v n -i(u n -i - v n -i), v n = v n (u n+ i - u n ). (43) 



Here, u n ~ v n ~ ^ . We computed a couple of symmetries for (43). One of them 
reads: 

G\ = u n _ x v n -\ + u n -\u n v n -\ + u n - 2 v n - 2 v n -i - v 2 _ 2 v n -i - 2u n _n^_ 1 
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G 2 



-u n v 2 l _ 1 + vl_i - u n u n+1 v n - u 2 n+1 v n + u n v 2 n + 2u n+ iv 2 n 
-v n - u n+2 v n v n+1 + v n vl +1 , 

U 2 V n - U 2 n+l V n + U n -iVn-lVn ~ V 2 n _ x V n - U n V 2 + U n+ \V 2 



-U n+ 2V n V n+ i + V n V n+1 . 

In [37], Suris investigated the integrability of 

tin = Vn [exp (y n+1 - y n ) - exp (y Tl 



Vn-l)] , 



(44) 



which is closely related to the classical Toda lattice (p2|). The same change of 

as 



variables as for (42) allows one to rewrite 



ii n = u n (v n - u„_i), i) n = v n (u n+ i - u n ). (45) 
Again, u n ~ v n ~ and we have computed three symmetries. Two of them are: 



4 2) 



n^U n + -UnUrt+iUn + U n V 2 - U n -\U n V n -\ - U 2 V n -i - U n V 2 n __ x , 
u n+l v n + Un+lV 2 + U n+ iV n V n+ i - U 2 n V n - U n V n -\V n - U n V 2 , 

u n v n + u 2 n u n+1 v n + u n u 2 n+1 v n + 2u n vl + 2u n u n+1 vl + u n v n 

+U n U n+ lV n V n+ l - v? n _ x u n v n -\ - U n -lU 2 n V n -l - U n V n -l 

-u n -iu n v n - 2 v n -i - 2u n -\u n v 2 n _ x - 2u^-u^_ 1 - u n v n _ 1 , 

u n+l v n - U n V n - U n ^iU n V n ^iV n - 2u 2 n V n -iV n - U n t^_ lUn - 2u^^ 

+2u n+1 v 2 - u n v n -ivl - u n v n + u n+1 v n + 2u 2 n+1 v n v n+ i 

+U n+ iU n+2 V n V n+ l + U n+ lVnV n+ l + U n+ lV n V 2 + 1 . 



5 The Integrability Package 

We now briefly describe the use of our Integrability Package, which has (among 
other things) the code for the computation of symmetries based on the algorithms 



in Sections 2 and 3. The Integrability Package is written in Mathematica [3S] 



syntax. Users are assumed to have access to Mathematica 3.0. All the necessary 



files are available in MathSource [15] including on-line help, documentation, and 
built-in examples. The corresponding files should be put in the appropriate places 
on your platform. Detailed instructions are given in the documentation. 

After launching Mathematica, type 
In[l]:= «Integrability' 



to read in the code. Doing so, you will get the following statement: 
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Loading init.m for Integrability from AddOns. 

For the purpose of symmetry computations, the functions PDESymmetries and 
DDESymmetries are available. 

Working with (|3^) as an example, the first two lines define the system (r = q*), 
whereas the third line will produce the three symmetries listed in (|33|): 

In[2]:= pdel:= D [q[x,t] ,t] - D[q[x,t] ,{x,2}] + 
2*q[x,t] ~2*r[x,t] == 0; 

In [3]:= pde2:= D[r[x,t],t] + D[r[x,t] ,{x,2}] - 
2*r [x,t] ~2*q[x,t] == 0; 

In [4] := PDESymmetries[{pdel,pde2},{q,r},{x,t},{4,6}, 
WeightRules->{Weight [q] ->Weight [r] }] 

Help about the functions and their options can be obtained by typing 

In [5]:= ??DDESymmetries 

DDESymmetries [eqn, u, {n, t}, R, opts] finds the symmetry with rank 
R of a differential-difference equation for the function u. 
DDESymmetries [{eqnl , eqn2, ...}, {ul, u2, ...}, {n, t}, R, opts] 
finds the symmetry of a system of differential-difference 
equations, where R is the rank of the first equation in the 
desired symmetry. DDESymmetries [{eqnl , eqn2, ...}, 
{ul, u2, ...}, {n, t}, {Rmax}, opts] finds the symmetries with 
rank through Rmax. DDESymmetries [{eqnl , eqn2, ...}, 
{ul, u2, ...}, {n, t}, {Rmin, Rmax}, opts] finds the symmetries 
with rank Rmin through Rmax. n is understood as the discrete 
space variable and t as the time variable. 

Attributes [DDESymmetries] = {Protected, ReadProtected} 

Options [DDESymmetries] = 

{WeightedParameters -> {}, WeightRules -> {}, 
MaxExplicitDependency -> 0, UndeterminedCoef f icients -> C} 

and 

In [6] := ??WeightedParameters 

WeightedParameters is an option that determines the parameters with 
weight. If WeightedParameters -> {pi, p2, ...}, 
then pi, p2, .... are considered as constant parameters with 
weight. The default is WeightedParameters -> {}. 
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Attributes [WeightedParameters] = {Protected} 

The option WeightedParameters is useful when working with systems that lack 
uniformity in rank. In such cases, the code tries to resolve the problem of lack of 
uniformity, and prints appropriate messages. If the code can not automatically 
resolve the problem it suggests the use of the WeightedParameters option. 
Therefore, one should not use the option WeightedParameters unless it is ex- 
plicitly suggested. For further descriptions of the functions and their options we 
refer to the documentation in Jl5f |. 



6 Other Software Packages 

Higher-order symmetries can be computed with prolongation methods and nu- 
merous software packages are available that can aid in the tedious computations 
inherent to such methods. An extensive review of software for Lie symmetry 
computations, including generalized symmetries, can be found in pl| 19|. 



With prolongation methods one generates and subsequently reduces and solves 
a determining system of linear homogeneous partial differential equations for the 
unknown higher-order symmetry. In many cases, due to the length and complexity 
of that system, the general solution is out of reach and one resorts to making a 
polynomial ansatz for the symmetry. 

Although restricted to polynomial higher-order symmetries, we believe that 
the method presented in this paper is much more straightforward. Furthermore, 
it does not require the application of prolongation methods or Lie algebraic tech- 
niques. 



To avoid retreading the surveys [18, 19], here, we restrict our discussion to 



symbolic packages that allow one to compute generalized symmetries of PDEs, 
as they were defined in Section 2.1. We are not aware of software for DDEs to 
calculate the type of symmetries defined in Section 3.1. 

Based on the alternative strategy discussed in Remark (ii) in Section 2, Ito's 



programs in REDUCE [22, 23] compute polynomial higher-order symmetries for 
systems of evolution equations that are uniform in rank (no weighted parameters 
can be introduced). Ito programs can not be used to computed symmetries that 
explicitly depend on x and t. 



In [10], Fuchssteiner et al. present an algorithm to compute higher-order 
symmetries of evolution equations. Their algorithm is based on Lie algebraic 
techniques and uses commutator algebra on the Lie algebra of vector fields. Their 
approach is different from the usual prolongation method in that no determining 
equations are solved. Instead, all necessary generators of the finitely generated 
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Virasoro algebra are computed from one given element by direct Lie algebra meth- 
ods. Their code is available in MuPAD. 

The REDUCE program FS for "formal symmetries" was written by Gerdt and 
Zharkov 12]. The code FS can be applied to polynomial nonlinear PDEs of evolu- 
tion type, which are linear with respect to the highest-order spatial derivatives and 
with non-degenerated, diagonal coefficient matrix for the highest derivatives. The 
algorithm in FS requires that the evolution equation are of order two or higher in 
the spatial variable. However, this approach does not require that the evolution 
equations are uniform in rank. With FS one cannot compute symmetries that 
depend explicitly on x and t. 

The PC package DELiA, written in Turbo PASCAL by Bocharov Q and 
co-workers, is a commercial computer algebra system for investigating differential 
equations using Lie's approach. The program deals with higher-order symmetries, 
conservation laws, integrability and equivalence problems. It has a special routine 
for systems of evolution equations. The program requires the presence of second 
or higher-order spatial derivative terms in all equations. 

Finally, Sanders and Wang [|33| have Maple and FORM software that aids in 
the computation of recursion operators. 



7 Conclusions 

We have implemented direct algorithms (in Mathematica) that allow the user to 
compute polynomial higher-order symmetries of polynomial systems of evolution 
and lattice equations. 

These algorithms are based on the dilation invariance of the given equations. 
Only minor modifications of our strategy lead to direct algorithms for conserved 
densities for systems of nonlinear PDEs [ 13 1 and nonlinear DDEs [[14], |i~6| ]. 



For systems that arise from a variational principle, conservation laws follow 
from higher-order symmetries (Noether's theorem) and vice versa. Currently, in 
our algorithms we are not exploiting such connections. 

In the future we will investigate generalizations of our methods to PDEs and 
DDEs in multiple space dimensions. We will also study the potential use of 
Lie-point symmetries other than dilation symmetries. Moreover, most recursion 
operators, which connect generalized symmetries, are also dilation invariant. We 
are extending our algorithms to the symbolic computation of recursion operators. 
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